/*  slmhhmul.cpp  by K.Tsuru  */
// function ID = 218    DRADIX, BRADIX
/**************************************************************************
SLong and SInteger classes
result = (huge m)*(huge n)  m >= n
Normal method using distribution law.
Using the FFT multiplication in which f*f=2f figures can be treated at the
maximum it provides the product with over f figures.
[Algorithm]
m.FFTSize() == m.FFTMaxArraySize();
It lets f = m.FFTMaxArraySize()/2 and takes F = R^f as the new radix, where
R=DRADIX or BRADIX.
mf, nf:the number of figures(blocks) in radix F.
m and n can be expressed as follows
    m = u(mf-1)*F^(mf-1)+u(mf-2)*F^(mf-2)+...+u(0)
    n = v(mf-1)*F^(nf-1)+v(nf-2)*F^(nf-2)+...+v(0)
Calculating u(i)*v(j) by FFT it is put in the (i+j)-th block of r.
When &m == &n(result = m*m) the calculation u(i)*v(j)(i > j) needs once.
**************************************************************************/
#ifndef SN_H
#include "sn.h"
#endif

void SLong::NHHMult(const SLong& m, const SLong& n, SLong& result){
    uint f = m.FFTSize()/2, d, k;
    int  i, j, mf, nf;
    bool m_sq = (&m == &n); // m == n, result = m*m ?
//figures of m = m.aHead + 1
    mf = (m.aHead+1)/f;     if( (m.aHead+1) % f ) mf++;
    nf = (n.aHead+1)/f;     if( (n.aHead+1) % f ) nf++;
/*
mf <= maxArraySize/f = 8
The "m.Size()" is a multiple of "f" because the "size" is a the power of two.
m >= n then maybe n.Size()<f.
If n.Size() <= f it directly use the value of "n" as "q".
*/
    uint qsz = (n.Size() <= f) ? 0 : f;
    if( (&result == &m) || (&result == &n)){
        m.SetError(m.SYNTAX_ERR, "NHHMult", 218);
    }
    result.valloc(m.aHead+n.aHead+2, 0);    //It allocates memory and initialize.

//worc area, The size of "t" is the same as that of "result".  p = q = t = 0
    SLong p(m.Type(), f), q(m.Type(), qsz), t(m.Type(), result.Size());
    const fType* mv = m.ReadFigures();
    const fType* nv = n.ReadFigures();
    fType* pv = p.figure.Elements();
    fType* qv = NULL;
    if(qsz) qv = q.figure.Elements();

#ifndef NDEBUG
m.figure(mf*f-1u);
#endif
    for(i = 0; i < mf; i++){
            p.SetSign(1);
       //It copys f figures to p down from the upper part of m.
            d = i*f;
            if(d > m.aHead) continue;      //   p = 0
            for(k = d ; k < d+f; k++) pv[k - d] = mv[k];
            p.CheckArray(218);  //It lets to know the LLMultFFT() the figure positions.
            if( !p.Sign(218) ) continue;
        for(j = 0; j < nf; j++){
            if(m_sq && (i < j)) continue;
            if(qsz){       // n.Size() > f
                d = j*f;
                if(d > n.aHead) continue;   // q = 0
                if( m_sq && (i == j) ){
                    t = LLMult(p, p);
                } else {
                    q.SetSign(1);
                    for(k = d ; k < d+f; k++) qv[k - d] = nv[k];
                    q.CheckArray(218);
                    if(q.Sign(218) == 0) continue;
                    t = LLMult(p, q);
                }
            } else {// nf = 1, n.Size() <= f, q = n
                t = LLMult(p, n);
            }
            if(m_sq && (i > j)) t = LsMult(t, 2);   // if i == j, t = p*p
            if(i+j) t.ShiftArray( (i+j)*f );        // t *= F^(i+j)
            result = LLAdd(result, t);              // result += t;
        }
    }
//A surplus of one figure was taken.It is a rare case.
    if( 2u*(result.aHead+1) <= result.figure.size() ) result.DoCutDown();
}